### Set Directory
try(setwd(dirname(rstudioapi::getActiveDocumentContext()$path)))
rm(list=ls())

### Load Package
library(stargazer)
library(dplyr)
library(tidyr)

### Load Data
data = read.csv('replication dataset.csv')

# Restrict sample to potential targets
# that have economic sanctions experienced
# threat and/or imposition during 1981-2005
data_main = data %>%
  filter(year %in% 1981:2005) %>%
  arrange(ccode,year) %>%
  group_by(ccode) %>%
  mutate(ever = ifelse(mean(threat)>0|mean(imposition)>0,1,0)) %>%
  ungroup() %>%
  filter(ever==1)

### First Stage OLS
m1 = lm(v2juncind.rev~v2juncind.rev.lag+hrdummy+nonhrdummy+polconiii+yr_sch_sec+open.log, data=data_main)
m2 = lm(v2x_corr~v2x_corr.lag+hrdummy+nonhrdummy+polconiii+yr_sch_sec+open.log, data=data_main)
m3 = lm(v2mecenefm.rev~v2mecenefm.rev.lag+hrdummy+nonhrdummy+polconiii+yr_sch_sec+open.log, data=data_main)

# Table 1
stargazer(m1,m2,m3,type='text')

### Second Stage Bootstrap

# Function to get fitted value
hat = function(d,m){
  X = d %>%
    select(names(coef(m))[-1])
  X = cbind(1,X)
  fit = as.matrix(X)%*%coef(m)
  return(fit)
}

# leave out observations with any missing data
d = data_main %>%
  select(v2juncind.rev.lag,v2x_corr.lag,v2mecenefm.rev.lag,
         hrdummy,nonhrdummy,polconiii,yr_sch_sec,open.log,
         tphysintrev,physintrev,
         gdppclog,polity2,civilwar,interwar) %>%
  drop_na()

# Function to bootstrap second stage
boot = function(x){
  d1 = d[sample(1:nrow(d),replace=T),]
  df1 =d1 %>%
    mutate(fit1 = hat(d1,m1)) %>%
    mutate(fit2 = hat(d1,m2)) %>%
    mutate(fit3 = hat(d1,m3))
  
  mb1 = lm(tphysintrev~physintrev+fit1+fit2+fit3
           +gdppclog+polity2+civilwar+interwar,
           data=df1)
  
  beta1 = coef(mb1)
  
  return(beta1)
}

# Bootstrapping (N=1000 iterations)
set.seed(2020)
N=1000
betas = replicate(N,boot())

# 95% Confidence Interval
result = t(apply(betas,1,function(x){
  coef = mean(x)
  lwr = sort(x)[floor(N*.025)]
  upr = sort(x)[floor(N*.975)]
  lwr2 = sort(x)[floor(N*.05)]
  upr2 = sort(x)[floor(N*.95)]
  return(cbind(coef,lwr,upr,lwr2,upr2))
}))

result = as.data.frame(result)

colnames(result) = c('coef','lwr','upr','lwr2','upr2')
result$id = 1:nrow(result)

# Table 2
round(result[,1:3],3)
nrow(d)

# Figure 2
par(mar=c(5,2,2,2))
plot(id~coef,data=result[3:5,],ylim=c(6,2),xlim=c(min(result[3:5,-6]),max(result[3:5,-6])),yaxt='n',xlab='Coefficient',ylab='',pch=c(15,16,17))
abline(v=0,lty='dashed')
for(i in 3:5){
  lines(y=c(i,i),x=c(result[i,2],result[i,3]))
}
text(y=result[3:5,6]-0.2,x=result[3:5,1],labels=c('Lower Court Dependence','Political Corruption','Government Censorship'),pos=4)

### Summary Statistics

# Table OA1
stargazer(as.data.frame(data_main %>%
                          select(physint,
                                 v2juncind.rev,
                                 v2x_corr,
                                 v2mecenefm.rev,
                                 hrdummy,
                                 nonhrdummy,
                                 polconiii,
                                 yr_sch_sec,
                                 open.log,
                                 gdppclog,
                                 polity2,
                                 civilwar,
                                 interwar)),type='text')
